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SUMMARY 

The  application  of  linear  programming  to  complex 


problems  gives  rise  to  the  handling  of  large  matrix  prob¬ 


lems.  The  simplex  method  -hae  proveg  to  be  the  most  exped¬ 
ient  of  any  tried  but,  since  this  is  essentially  Gaussian 
elimination,  numerical  techniques  for  handling  it  have 
definite  limitations,  both  with  regard  to  running  time  and 
memory  size. 

RAND's  simplex  codes  for  the  IBM  701  hare  been 
designed,  through  revisions  to  the  method  and  special  com¬ 
piling  routinen,  to  maintain  great  precision  and  facilitate 
checking,  to  be  adaptable  to  many  variations,  and  to  be  as 
automatic  in  operation  as  possible.  Several  quite  large 
models  of  considerable  variety  are  discussed  along  with 
the  prospects  for  future  improvements,  both  in  theory  and 
in  coding  techniques.  (  )  ^ 


! 
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COMPUTINQ  EXPERIENCE  WITH  LINEAR  PROORAMMINQ 

AND  ITS  VARIANTS 


A  great  deal  has  been  said  and  written  in  the  past  two 
or  th2*ee  years  about  the  technique  known  as  linear  program¬ 
ming  which  uses  a  mathematical  model  to  represent  a  real 
situation.  When  one  undertakes  to  study  a  complex  real-life 
problem  by  such  a  model,  he  quickly  finds  that  the  paradise 
of  optimality  is  guarded  by  three  monsters  called  Conceptu¬ 
alization,  Porraulatlon,  and  Solution,  Any  apparent  defeat 
of  this  trio  turns  out  to  be  only  temporary  since  they  seem 
capable  of  escaping  any  conceivable  kind  of  shackles.  How¬ 
ever,  they  are  not  completely  Invulnerable  and  linear  pro¬ 
gramming  is  Interesting  and  usefiil  precisely  because  it  is 
linear  and  hence  amenable  to  relatively  simple  mathematical 
procedures.  Since  the  formal  definition  refers  to  the  model 
rather  than  to  its  antitype,  it  is  meaningful  only  after  one 
has  had  some  exposure  to  the  method.  It  will  be  assumed 
that  the  reader  is  familiar  with  the  form  of  such  a  model 
and  with  the  more  common  terms  of  the  subject  ,* 

The  old  cliche  about  not  seeing  the  forest  for  the  trees 
applies. especially  well  to  the  question  of  conceptualization. 
Oftentimes,  one  has  lived  to  close  to  a  tough  problem  for  no 
long  that  he  has  a  n«ntal  block  against  a  new  approach  which 
appears  to  make  ruthless  simplifications.  Yet  such  a  step 
is  often  necessary  in  applying  a  new  technique.  Clearly,  it 
will  sometimes  be  necessary  to  oversimplify,  or  at  least  to 


*  Bracketed  numbers  refer  to  references  at  end  of  paper. 
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take  a  different  view  of  some  factorsi  In  order  to  fit  a  prob« 
lem  Into  a  linear  context.  Furthermore,  one  must  usually  be 
very  selective  In  chocsing  the  conditions  that  one  Is  going 
to  Impose  on  a  model.  One  cannot  hope  to  take  account  of, 
say,  every  facet  of  the  dally  activities  of  an  Industry  In  a 
quantitative  statement.  Nevertheless,  In  spite  of  these  con¬ 
ceptual  difficulties,  the  range  of  problems  for  which  useful 
results  can  be  obtained  Is,  I  believe,  surprisingly  large. 

Once  one  has  conceived  of  a  way  to  represent  a  realistic 
number  of  co..stralnts  without  running  their  number  up  beyond 
computational  capabilities,  he  Is  faced  with  the  task  of 
actually  formulating  the  numerical  model.  This  process,  It 
should  be  noted,  forces  someone  to  acquire  an  Intimate  and 
detailed  knowledge  of  the  complex  being  studied.  Our 
customers  are  often  able  to  Interpret  the  results  of  numerical 
computations  very  lucidly,  often  while  the  problem  Is  mmnlng. 
This  sort  of  Insight  is  valuable  In  Itself  so  long  as  one  dis¬ 
tinguishes  between  the  model  and  the  real  world. 

Vfhen  the  numerical  model  Is  set  down,  there  Is  the  prob  - 
lem  of  how  to  handle  the  computations.  Furthermore,  one 
optimal  solution  Is  quite  likely  to  be  Insufficient  to  answer 
all  the  formulator's  questions.  (In  fact,  we  have  come  to 
view  with  suspicion  any  problem  which  Is  to  be  run  only  once.) 
If  this  Is  known  In  advance,  we  can  often  suggest  changes  In 
the  model,  or  at  least  In  the  way  of  looking  at  It,  which 
will  render  these  continuing  computations  easier  to  perform. 

For  example,  in  some  problems  which  were  not  too  large.  It 
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was  discovered  that  It  would  be  necessary  to  add  additional 
restraints  after  each  optimal  solution  was  computed.  By 
dualizing  the  problemj  It  was  possible  to  hold  the  number 
of  restraints  fixed  throughout  and  merely  Increase  the  number 
of  activity  vectors  from  time  to  time,  proceeding  from  one 
optimal  solution  to  the  next  without  starting  all  over  with 
a  larger  system. 

It  should  be  explained,  Incidentally,  that  all  the 
problems  whose  solutions  we  have  been  concerned  with  have 
been  solved  with  the  simplex  method  or  some  variant  thereof. 
The  number  of  these  variants  Is  constantly  growing  and  It 
appears  that  the  trend  will  be  to  specialized  procedures  for 
problems  classified  by  types,  at  least  for  large  problems 
[^,5]»  However,  the  present  discussion  pertains  to  general 
methods. 

The  first  problem  for  which  extensive  runs  were  made  at 
RAND  was  a  petroleum  blending  study  [6,  7]  •  The  calculations 
were  carried  out  on  an  IBM  CPC  Model  II  with  a  special  set-up. 
The  product  form  for  the  inverse  of  the  basis  matrix  was 
first  used  In  this  set-up  to  speed  up  operation  [8,  9] •  This 
form  has  been  carried  over  into  all  our  later  simplex  codes 
for  the  IBM  701  and  has  proved  to  be  an  extremely  versatile 
tool  for  working  with  matrices  and  their  Inverses  on  high¬ 
speed  conputers . 

We  also  performed  what  has  become  known  as  parametric 
programming,  or  PLP,  on  this  first  problem.  PLP  can  be 
divided  Into  three  major  types, according  as 
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(1)  the  elements  of  the  right-hand  side 

(2)  the  coefficients  of  the  objective  function,  or 

(3)  elements  of  the  restraint  matrix 

are  expressed  as  linear  functions  of  a  single  parameter 
Chap.  4].  All  three  of  these  were  applied,  although 
quite  crudely  at  the  time.  Parametrizing  elements  of  the 
matrix  Is  not  a  very  satisfactory  device  since  one  must 
simultaneously  maintain  optimality  and  feasibility  and  also 
avoid  passing  through  a  point  of  singularity.  As  a  result, 
the  calculations  become  so  Involved  that  the  cost  of  per¬ 
forming  them  eats  up  the  desired  saving.  However,  alternate 
techniques  for  making  changes  within  the  matrix  have  been 
devised,  mainly  for  correcting  errors  found  after  much  com¬ 
puting  time  has  been  invested.  It  Is  now  very  seldom  that 
a  problem  must  be  restarted  from  the  beginning. 

We  have  been  able  to  develop  the  parametrlzatlon  of  the 
right-hand  side  Into  a  fully  automatic  procedure,  and  It  has 
been  used  extensively.  Variations  In  the  optimizing  form  are 
not  as  neat  computing-wise .  The  reason  for  this  Is  that  one 
must  repeatedly  alter  a  set  of  n  numbers,  where  n  Is  the 
number  of  variables,  Instead  of  m  numbers,  where  m  Is  the 
number  of  constraints  and  typically  much  smaller  than  n. 

This  Is  at  odds  with  the  basic  phlxosophy  of  our  code.  How¬ 
ever,  here  again  alternate  devices  have  been  perfected  which 
In  practice  accomplish  virtually  the  same  result.  Our  latest 
code,  for  example,  will  optimize  successively  and  automatic¬ 
ally  on  a  multiplicity  of  objective  functions.  These  forms 
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can  be  made  to  differ  by  non-uniform  ratios  so  that  in  a 
sense  It  is  more  general  than  actual  PLP.  The  optimal  solu¬ 
tion  for  the  first  form  is  used  as  a  starting  point  for  the 
second,  the  optimal  solution  for  the  second  as  a  starting 
point  for  the  third,  and  so  on.  On  the  other  hand,  this 
method  does  not  permit  one  to  solve  for  the  critical  values 
of  a  single  parameter  automatically.  We  can  obtain  these  one 
at  a  time  without  too  much  effort  but  they  are  not  usually 
required. 

It  was  believed  that  an  extremely  fast  code  using  fixed 
point  arithmetic  could  be  programmed  for  the  701  which  would 
have  been  capable  of  handling  systems  with  100  restraint  equa¬ 
tions.  Such  a  program  was  Indeed  coded  but  it  proved  to  be 
impractical  due  to  the  fast  accumulation  of  round-off  error. 

A  second  version  was  patched  up  from  the  first  but  it  gave  no 
better  results.  \  j  then  decided  to  go  to  double-precision 
floating-point  arithmetic  for  the  third  code  and  this  did  per¬ 
mit  us  to  heuidle  100  order  systems  successfully.  Since  the 
simplex  method  is  essentially  selective  Qausslan  elimination, 
it  is  mandatory  that  one  maintain  extreme  precision  even 
though  the  input  data  may  not  be  very  accurate  Itself.  An 
error  analysis  is  small  comfort  when  one  must  require  the 
machine  to  make  decisions  on  the  basis  of  intermediate  results. 
Thus  we  feel  that  the  extra  computing  time  is  a  necessary  cost 
for  a  successful  code. 

However,  as  the  code  became  more  elaborate,  it  became 
almost  impossibly  Interwoven  so  that  a  fourth  code  was  devised 
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whlch  emphasized  clean  logical  structure  [lo]  .  The  present 

code,  which  is  the  fifth,  is  only  a  slight  refinement  of  the 

fourth, emphasizing  simplicity  of  operation.  Space  has  been 

provided  for  up  to  25O  restraint  equations  in  almost  any 

reasonable  number  of  variables  --  1553  being  the  largest  to 

\ 

date  --  but  solution  times  are  somewhat  high  for  systems  with 
more  than  15O  equations.  However,  the  various  parts  of  the 
code  are  designed  so  that  they  can  be  coupled  together  to 
form  new  routines  almost  at  will.  This  is  accomplished  by 
using  a  compiling  routine  which  puts  tcj^ether  subroutines 
Imroedlately  prior  to  their  use  while  a  problem  is  running.  A 
master  control  region  activates  the  compiler  througli  a  combina¬ 
tion  link  and  pseudo^instructlon.  Naturally,  we  have  auxili¬ 
ary  routines  for  use  with  the  system  to  facilitate  re-start- 
Ing,  picking  up  after  a  machine  failure,  or  altering  the 
model,  as  well  as  to  make  certain  side  computations  sometimes 
required.  In  addition,  a  composite  algorithm  built  Into  the 
basic  routine  takes  over  automatically  If  any  Infeaslbllltlfss 
exist  In  the  Initial  solution.  This  Is  often  used  In  place 
of  the  dual  simplex  algorithm  [ll]  since  the  latter  is  slower 
in  operation  and  assumes  optimality  as  a  starting  point.  In 
short,  we  have  attempted  to  provide  for  a  variety  of  problems 
In  two  ways:  first,  by  having  a  basic  routine  which  is  as 
general  and  elaborate  as  practicable;  and  second,  by  allowing 
a  variety  of  programs  to  be  constructable  in  our  coding  sys¬ 
tem  with  a  mlnlmuir.  of  effort,  sort  of  a  programmer's  set  of 
"pluggable  units". 
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The  initial  petroleum  model  was  followed  by  three 
others.  The  first  was  a  105  order  system  which  represented 
in  considerable  detail  the  U.  S.  refining  and  crude  oil  pro¬ 
ducing  industry  [l2]  .  This  model,  like  the  first,  was  formu¬ 
lated  by  Dr.  Alan  S.  Manne.  After  obtaining  an  initial 
optimal  solution,  the  right-hand  side  was  parametrized  in  two 
ways  and  two  broken-line  trade-off  curves  were  obtained  be¬ 
tween  jet  and  other  refinery  products.  Using  one  para¬ 

meter,  it  was  possible  to  trace  through  the  effects  of  vary¬ 
ing  the  jet  fuel  requirement.  The  other  parameter  varied 
the  availability  of  certain  refinery  equipment  and  was  used 
only  to  get  from  one  curve  to  the  other.  Each  turning  point 
on  each  curve  represented  the  solution  of  a  linear  program¬ 
ming  problem.  These  critical  values  of  the  parameter  are 
determined  automatically  in  this  type  of  PLP  and  are  some¬ 
times  of  considerable  Interest,  as  in  this  case. 

This  model  was  aggregated  to  cut  down  the  order  of  the 
system  and  parallel  computations  made  to  determine  the  sensi¬ 
tivity  to  detail.  With  appropriate  adjustments,  the  aggre¬ 
gated  version  was  then  blown  up  into  an  inter-regional  model 
involving  transportation  as  well  as  production  facilities. 
This  model  has  195  restraint  equations  in  more  than  1500 
variables.  The  computations  are  proceeding  well  and  some  PLP 
is  planned.  In  handling  this  many  variables,  however,  it  is 
wise  to  progressively  sub-optimlze  on  nested  subsets  of  the 
variables,  a  procedure  which  has  been  followed  on  several  oc¬ 


casions  . 
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Thls  last  technique  Is  well  illustrated  by  a  large  air¬ 
transport  scheduling  model  also  developed  by  Dr.  Mannt  [l3]  • 
The  system  was  so  large  that  the  generation  of  the  matrix 
Itself  was  done  by  machine.  Three  different  cases  were  run, 
all  with  94  restraint  equations  plus  a  minimizing  form,  and 
with  the  following  numbers  of  variables: 

Case  1,  1476  variables. 

Case  2,  1494  variables  Including  those  of  case  1, 

Case  3#  1553  variables  Including  all  but  one  of  those 
of  case  2. 

First,  however,  a  set  of  94  vectors  which  were  known  to  be 
linearly  independent  were  pre-selected  to  form  an  Initial 
basis  which  was  Inverted  In  73  Iterations,  there  being  21 
unit  vectors  In  the  set.  The  Inversion  Iterations  go  much 
faster  than  normal  simplex  cycles,  and  this  Inversion  took 
about  30  minutes.  This  gave  a  solution  which  was  Infeasible 
In  8  variables.  However,  by  using  an  expanded  set  of  123 
vectors.  Including  the  original  94,  the  composite  algorithm 
reduced  this  to  one  small  Infeasiblllty  In  27  Iterations  of 
which  3  were  repetitive  ones  to  correct  some  errors  which  had 
been  detected  In  the  data. 

The  other  1353  vectors  for  case  1  were  then  made  avail¬ 
able  and  a  feasible  solution  obtained  In  8  cycles,  meanwhile 
Improving  the  minimizing  form.  After  75  more  cycles.  It  be¬ 
came  advisable  to  re-invert  the  current  basis  to  reduce  the 
length  of  the  transformations  tape.  This  required  89  of  the 
fast  Iterations,  and  after  57  more  simplex  cycles,  optimality 
was  obtained  for  case  1.  Cases  2  and  3  were  picked  up  In  sue 
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ceaslon  In  a  olrallar  manner.  The  total  number  of  Iterations 
for  all  cases  was  739*  of  which  359  were  fast  Inversion  cycles. 
Thus  380  Iterations  were  performed  to  obtain  new  results  In 
solving  all  three  problems,  quite  good,  I  believe,  consider¬ 
ing  the  very  large  number  of  variables. 

Harry  Markowitz  has  used  the  simplex  codes  at  RAND  to 
solve  problems  In  which  the  variables  are  allowed  to  take  on 
only  Integral  values,  or  what  might  be  called  Diophantlne 
linear  programming  Cl^] •  This  is  not  done  by  a  single  run, 
but  by  adding  additional  constraints  from  time  to  time  on  the 
basis  of  previous  optimal  solutions.  A  process  of  elimina¬ 
tion  Is  used  to  reduce  the  number  of  admissible  variables, 
not  by  omitting  them,  but  by  forcing  them  above  or  below  cer¬ 
tain  bounds  with  the  added  constraints.  The  bounds  are  com¬ 
puted  on  the  side. 

The  code  has  been  used  for  fitting  t  observations  to 
m  variables,  minimizing  the  sum  of  the  absolute  values  of  the 
errors.  This  required  using  2t  +  2m  variables  as  It  was  set 
up.  This  process  could  be  made  more  efficient  by  a  special 
routine  but  we  have  not  taken  time  yet  to  code  It.  It  would 
differ  from  the  regular  simplex  method  mainly  In  the  way  in 
which  the  variable  to  be  eliminated  Is  chosen  on  each  cycle. 

As  Is  fairly  well  known,  the  simplex  method  Is  a  good 
way  to  solve  linear  systems  since  It  often  gives  valuable  In¬ 
formation  even  when  no  solution  exists.  We  have  also  found 
the  code  very  good  for  Inverting  badly  conditioned  matrices. 

If  the  order  Is  not  too  big  --  say  50,  although  the  practical 
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value  of  such  Inverses  Is  open  to  question. 

We  have  recently  been  working  on  a  nutrition  problem 
which  Is  unusual  In  that  there  are  only  l6  restraint  equa¬ 
tions,  3  of  these  somewhat  artificial,  but  several  hundred 
variables.  Furthermore,  the  matrix  Is  very  full.  In  fact  It 
contains  more  non-zero  entries  than  any  other  problem  we  have 
run.  It  Is  b  ing  used  to  study  objective  functions  and  other 
aspects  cf  model  formulation,  with  some  Interesting  and  bizarre 
diets  as  a  by-product. 

We  are  already  preparing  a  704  code  which  will  be  essen¬ 
tially  the  same  as  our  present  701  code.  However,  two  new 
features  are  planned.  In  the  701,  the  original  matrix  Is 
carried  In  condensed  form,  that  Is,  only  non-zero  elements 
are  recorded  and  these  are  Indexed.  We  plan  to  do  the  same 
with  the  transformations  on  the  704.  This,  together  with 
the  higher  speed  of  the  704  Itself,  should  enable  us  to  handle 
200  or  perhaps  250  order  systems  In  very  acceptable  times .  We 
now  feel  that  this  Is  €in  upper  limit  for  a  general  routine, 
not  only  because  of  time,  but  also  because  round-off  error 
again  becomes  troublesome  at  this  point  even  with  double  pre¬ 
cision  numbers. 

The  other  feature  is  variables  with  upper  bounds.  This 
can  be  handled  without  much  cost  either  complexity  or  time 
If  certain  conventions  are  adopted.  Since  several  restraint 
equations  are  sometimes  used  merely  for  bounding  certain 
activities,  this  will  allow  us  to  handle  effectively  bigger 
systems  for  such  problems. 
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We  have  a  code  for  RAND's  ovm  machine «  the  JOHNNIAC« 
which  uses  a  different  elimination  method,  developed  by 
Harry  Ilarkowltz,  for  inverting  a  matrix  [15] .  It  is  designed 
to  maintain  as  many  zeros  as  possible  in  the  representation 
of  the  Inverse  of  a  sparse  matrix,  as  linear  progrcunming  bases 
typically  are.  The  code  has  not  yet  been  in  operation  long 
enough  to  assess  its  value  but  the  method  is  compatible  with 
the  simplex  algorithm.  Beyond  this,  we  know  of  no  other  ways 
to  make  a  general  routine  handle  larger  systems  efficiently. 
Incidentally,  it  seems  odd  that  no  one,  so  far  as  I  know,  has 
before  attempted  to  solve  a  linear  system  on  a  machine  using 
the  same  rules  of  thumb  that  one  would  use  In  solving  by 
hand,  l.e.  to  eliminate  first  on  the  rows  and  columns  having 
the  fewest  elements. 

Alan  Manne  noted  that  his  air-transport  model  had  a 
matrix  structure  which  was  a  special  case  of  block  triangu¬ 
larity.  Upon  Investigation,  he  found  that  many  problems  fell 
into  this  pattern.  He  and  Dr.  Oeorge  Dantzlg  worked  out  a 
means  of  exploiting  this  structure  and  they,  Ted  Robacker, 
and  myself  are  now  working  out  the  details  of  a  machine 
method  to  take  advantage  of  magnetic  tape  characteristics. 
Manne  Is  planning  to  spend  considerable  time  Investigating 
such  semi-spec lallzed  methods  and  this  is  currently  the  most 
promising  means  of  handling  larger  systems. 
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